Investigating the modification of the DRHODZ variable in the MSIS.f90 subroutine.

DRHODZ is the part of the calculation of the partial derivative of mass density in the D matrix of the variational equations.

This notebook looks at how updating the calculation impacts the DRHODZ term.

Want to show:

  1. MSIS00

    1. original DRHODZ

    2. updated DRHODZ

  2. MSIS2

    1. original DRHODZ

    2. updated DRHODZ

Import base parameters

[1]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
from plotly.subplots import make_subplots
import plotly.express as px

import copy
import sys
import numpy as np
[2]:
%load_ext autoreload
%autoreload 2

import sys
sys.path.insert(0, '/data/geodyn_proj/pygeodyn/')
from pygeodyn_ReadOutput import pyGeodyn_Readers

run_params = {}
run_params['run_ID']           = 'RUN'
run_params['satellite']        = 'starlette'
run_params['den_model']        = None
run_params['empirical_accels'] = False
run_params['SpecialRun_name']  = None
run_params['options_in']       = {'DRHODZ_update':True}
run_params['verbose']          = False

run_params['arc'] = [ '030914_2wk',
                      '030928_2wk',
                      '031012_2wk',
                      '031026_2wk',
                      '031109_2wk',
                      '031123_2wk',
                      '031207_2wk'
                     ]


Load different versions

[3]:
#---MSIS2-DRHODZ-ON---------------------------------
params_msis2_DrhodzOn  = copy.deepcopy(run_params)
params_msis2_DrhodzOn['den_model'] = 'msis2'
params_msis2_DrhodzOn['SpecialRun_name']  = ''
params_msis2_DrhodzOn['options_in']       = {'DRHODZ_update':True}
Obj_msis2_DrhodzOn = pyGeodyn_Readers( params_msis2_DrhodzOn )
# Obj_msis2_DrhodzOn.getData_OnAllArcs()


#---MSIS2-DRHODZ-OFF---------------------------------
params_msis2_DrhodzOff  = copy.deepcopy(run_params)
params_msis2_DrhodzOff['den_model'] = 'msis2'
params_msis2_DrhodzOff['SpecialRun_name']  = '_drhodzOrig'
params_msis2_DrhodzOff['options_in']       = {'DRHODZ_update':False}
Obj_msis2_DrhodzOff = pyGeodyn_Readers( params_msis2_DrhodzOff )
# Obj_msis2_DrhodzOff.getData_OnAllArcs()


#---MSIS00-DRHODZ-ON---------------------------------
params_msis00_DrhodzOn  = copy.deepcopy(run_params)
params_msis00_DrhodzOn['den_model'] = 'msis00'
params_msis00_DrhodzOn['SpecialRun_name']  = ''
params_msis00_DrhodzOn['options_in']       = {'DRHODZ_update':True}
Obj_msis00_DrhodzOn = pyGeodyn_Readers( params_msis00_DrhodzOn )
# Obj_msis00_DrhodzOn.getData_OnAllArcs()


#---MSIS00-DRHODZ-OFF---------------------------------
params_msis00_DrhodzOff  = copy.deepcopy(run_params)
params_msis00_DrhodzOff['den_model'] = 'msis00'
params_msis00_DrhodzOff['SpecialRun_name']  = '_drhodzOrig'
params_msis00_DrhodzOff['options_in']       = {'DRHODZ_update':False}
Obj_msis00_DrhodzOff = pyGeodyn_Readers( params_msis00_DrhodzOff )
# Obj_msis00_DrhodzOff.getData_OnAllArcs()


Calling pygeodyn with multiple arcs...
Calling pygeodyn with multiple arcs...
Calling pygeodyn with multiple arcs...
Calling pygeodyn with multiple arcs...

Read Data into Objects

[4]:
Obj_msis2_DrhodzOn.getData_OnAllArcs()

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st030914_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st030928_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031012_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031026_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031109_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031123_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031207_2wk.goco05s
[5]:
Obj_msis2_DrhodzOff.getData_OnAllArcs()

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st030914_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st030928_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st031012_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st031026_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st031109_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st031123_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st031207_2wk.goco05s
[6]:
Obj_msis00_DrhodzOn.getData_OnAllArcs()

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st030914_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st030928_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031012_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031026_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031109_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031123_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031207_2wk.goco05s
[7]:
Obj_msis00_DrhodzOff.getData_OnAllArcs()

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st030914_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st030928_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st031012_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st031026_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st031109_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st031123_2wk.goco05s

     File path:
     Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st031207_2wk.goco05s

Plot Results

[8]:
# Obj_msis00_DrhodzOff.Density['']

# Obj_msis00_DrhodzOff.__dict__['Density']
# X = obj_m1.__dict__['Density'][arc]['Date'][::decimate]
# Y = obj_m1.__dict__['Density'][arc]['drhodz (kg/m**3/m)'][::decimate]

[9]:
%load_ext autoreload
%autoreload 2

import sys
sys.path.insert(0, '/data/geodyn_proj/pygeodyn/')
from pygeodyn_Analysis import plot_drhodz


fig = make_subplots(rows=2, cols=1   )

fig = plot_drhodz(fig, Obj_msis00_DrhodzOff,Obj_msis00_DrhodzOn, 0, 200, 'old drhodz','new drhodz' )
# fig = plot_drhodz(fig, , 1, 200, 'DRHODZ update On')

iplot(fig)
The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
[10]:
%load_ext autoreload
%autoreload 2

import sys
sys.path.insert(0, '/data/geodyn_proj/pygeodyn/')
from pygeodyn_Analysis import plot_drhodz


fig = make_subplots(rows=2, cols=1   )

fig = plot_drhodz(fig, Obj_msis2_DrhodzOff,Obj_msis2_DrhodzOn, 0, 200, 'old drhodz','new drhodz' )

iplot(fig)
The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
[ ]:

[ ]:

[ ]:

[11]:
# # params_msis00_origDrhoDz
# params_msis2_origDrhoDz
[12]:
# from pygeodyn_ReaderDevelop import pyGeodyn_Readers

# Obj_msis00_origDrhoDz = pyGeodyn_Readers( params_msis00_origDrhoDz )
# Data_msis00_origDrhoDz = Obj_msis00_origDrhoDz.getData_UserChoose(['Density'])

# # Obj_msis2_origDrhoDz = pyGeodyn_Readers( params_msis2_origDrhoDz )
# # Data_msis2_origDrhoDz = Obj_msis2_origDrhoDz.getData_UserChoose(['Density'])

# Obj_msis00_updateDrhoDz = pyGeodyn_Readers( params_msis00_updateDrhoDz )
# Data_msis00_updateDrhoDz = Obj_msis00_updateDrhoDz.getData_UserChoose(['Density'])

# # Obj_msis2_updateDrhoDz = pyGeodyn_Readers( params_msis2_updateDrhoDz )
# # Data_msis2_updateDrhoDz = Obj_msis2_updateDrhoDz.getData_UserChoose(['Density'])

[13]:




# import plotly.graph_objects as go
# import numpy as np
# from plotly.offline import plot, iplot
# # %matplotlib inline
# from plotly.subplots import make_subplots
# import plotly.express as px
# # import plotly
# # plotly.offline.init_notebook_mode(connected=True)

# col1 = px.colors.qualitative.Plotly[0]
# col2 = px.colors.qualitative.Plotly[1]
# col3 = px.colors.qualitative.Plotly[2]
# col4 = px.colors.qualitative.Plotly[3]

# data_nums_1 = 1000
# data_nums_2 = 100

# fig = make_subplots(rows=2, cols=1,
#     subplot_titles=("Direct Comparison","Percent Difference in Original and Edited DRHODZ"),
#                    )

# Obj00 = Data_msis00_origDrhoDz.Density
# Obj01 = Data_msis00_updateDrhoDz.Density
# # Obj10 = Data_msis2_origDrhoDz.Density
# # Obj11 = Data_msis2_updateDrhoDz.Density


# # for arcs in arc_list:
# fig.add_trace(go.Scatter(x=Obj00['Date'][::data_nums_2],
#                          y=np.abs(Obj00['drhodz (kg/m**3/m)'][::data_nums_2]),
# #                             name=name_m1,
#                              mode='markers',
#                             marker=dict(color=col1,
#                             size=4,
#                             ),
#                            ),
#                             row=1, col=1,
#                            )


# fig.add_trace(go.Scatter(x=Obj01['Date'][::data_nums_2],
#                          y=np.abs(Obj01['drhodz (kg/m**3/m)'][::data_nums_2]),
# #                             name=name_m1,
#                              mode='markers',
#                             marker=dict(color=col2,
#                             size=4,
#                             ),
#                            ),
#                             row=1, col=1,
#                            )


# # fig.add_trace(go.Scatter(x=Obj10['Date'][::data_nums_2],
# #                          y=np.abs(Obj10['drhodz (kg/m**3/m)'][::data_nums_2]),
# # #                             name=name_m1,
# #                              mode='markers',
# #                             marker=dict(color=col3,
# #                             size=4,
# #                             ),
# #                            ),
# #                             row=1, col=1,
# #                            )


# # fig.add_trace(go.Scatter(x=Obj11['Date'][::data_nums_2],
# #                          y=np.abs(Obj11['drhodz (kg/m**3/m)'][::data_nums_2]),
# # #                             name=name_m1,
# #                              mode='markers',
# #                             marker=dict(color=col4,
# #                             size=4,
# #                             ),
# #                            ),
# #                             row=1, col=1,
# #                            )

# A0 = (Obj00['drhodz (kg/m**3/m)'][::data_nums_2])
# B0 = (Obj01['drhodz (kg/m**3/m)'][::data_nums_2])
# C0 = ((B0-A0)/A0)*100

# fig.add_trace(go.Scatter(x=Obj00['Date'][::data_nums_2],
#                          y=C0,
#                          name = '% diff',
#                          mode='markers',
#                             marker=dict(color='green',
#                             size=5,
#                             ),
# #                           showlegend=False,
#                            ),
#                            row=2, col=1,
#                            )
# # A1 = (Obj10['drhodz (kg/m**3/m)'][::data_nums_2])
# # B1 = (Obj11['drhodz (kg/m**3/m)'][::data_nums_2])
# # C1 = ((B-A)/A)*100

# # fig.add_trace(go.Scatter(x=Obj10['Date'][::data_nums_2],
# #                          y=C,
# #                          name = '% diff',
# #                          mode='markers',
# #                             marker=dict(color='red',
# #                             size=3,
# #                             ),
# # #                           showlegend=False,
# #                            ),
# #                            row=2, col=1,
# #                            )



# fig.update_layout(
#     title="Modifying DRHODZ-- Comparison",
#     )
# fig.update_yaxes(type="log", exponentformat= 'power', row=1, col=1)
# # fig.update_yaxes( exponentformat= 'power', row=2 , col=1)

# fig.update_yaxes(   exponentformat= 'power',)
# fig.update_layout(legend= {'itemsizing': 'constant'})

# fig.update_yaxes(title_text="kg/m^4", row=1, col=1)
# fig.update_yaxes(title_text="%", row=2, col=1)

# fig.update_xaxes(title_text="Date", row=1, col=1)
# fig.update_xaxes(title_text="Date", row=2, col=1)
# fig.update_layout(legend= {'itemsizing': 'constant'})
# #     fig.update_xaxes(tickangle=45)
# fig.update_layout(
#     autosize=False,
#     width=900,
#     height=700,
# )
# fig.update_layout(
#     font=dict(
# #             family="Courier New, monospace",
#         size=14,
#              ),)


# iplot(fig)

[14]:
# import os

# # if not os.path.exists("msis_update_images"):
# #     os.mkdir("msis_update_images")

# fig.write_image("plot_drhodz_msis00_eval.png")

[ ]:

[ ]:

[ ]:

[ ]: